Measurements of the Yield Stress in Prictionless Granular Systems 
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We perform extensive molecular dynamics simulations of 2D frictionless granular materials to 
determine whether these systems can be characterized by a single static yield shear stress. We 
consider boundary-driven planar shear at constant volume and either constant shear force or constant 
shear velocity. Under steady flow conditions, these two ensembles give similar results for the average 
shear stress versus shear velocity. However, near jamming it is possible that the shear stress required 
to initiate shear flow can differ substantially from the shear stress required to maintain flow. We 
perform several measurements of the shear stress near the initiation and cessation of flow. At fixed 
shear velocity, we measure the average shear stress Ej;„ in the limit of zero shear velocity. At fixed 
shear force, we measure the minimum shear stress Ey/ required to maintain steady flow at long times. 
We find that in finite-size systems E^/ > Ey„, which implies that there is a jump discontinuity in 
the shear velocity from zero to a finite value when these systems begin flowing at constant shear 
force. However, our simulations show that the difference Ey/ — Eyi,, and thus the discontinuity in 
the shear velocity, tend to zero in the infinite system size limit. Thus, our results indicate that in 
the large system limit, frictionless granular systems are characterized by a single static yield shear 
stress. We also monitor the short-time response of these systems to applied shear and show that 
the packing fraction of the system and shape of the velocity profile can strongly influence whether 
or not the shear stress at short times overshoots the long-time average value. 

PACS numbers: 47.50.-|-d 83.10.Mj 83.50.-v 45.70.Mg, 



I. INTRODUCTION 

The static yield shear stress, or similarly the static 
shear modulus, is an important material property that 
distinguishes soHds from liquids 0. Solids possess a 
nonzero static yield shear stress, while it vanishes for liq- 
uids. Solids are able to resist applied shear stresses below 
the yield shear stress, but plastic flow occurs when shear 
stresses larger than the yield shear stress are applied. In 
contrast, liquids flow when any finite shear stress is ap- 
plied. 

Disordered materials such as molecular and colloidal 
glasses, static granular materials, and concentrated emul- 
sions also possess a nonzero yield shear stress. However, 
it is difficult to determine precisely the yield shear stress 
in these amorphous systems since they often display non- 
linear and spatially nonuniform response, for example 
creep fiow, intermittent dynamics, and shear localization, 
when shear stress is applied. The value of the yield stress 
in these amorphous systems can also depend on how it 
is measured. For example, the yield shear stress required 
to generate steady flow in an originally unsheared sys- 
tem may differ significantly from a measure of the yield 
stress obtained by approaching the static state by slowly 
decreasing the shearing velocity. The yield shear stress 
may also depend strongly on how the system was pre- 
pared. For example, it has been shown that the yield 
shear stress is sensitive to the age and strain history in 
glassy systems 0] and the construction history 0, 0| and 
micro-structural details [5| in granular materials. 



There have been a number of recent computational in- 
vestigations of the transition from static to flowing states 
in granular and glassy systems. For example, measure- 
ments of the yield shear stress or static shear modulus 
have been conducted as a function of packing fraction in 
model foams 0, emulsions 01 , and frictionless granular 
materials and as a function of temperature and strain 
rate in dense Lennard- Jones glasses |^ liol ITU , metal- 
lic glasses l^, and polymer glasses |l3l Il4j. However, 
an important question that has not been adequately ad- 
dressed by these previous studies is whether or not there 
is a unique measure of the static yield shear stress in 
amorphous granular and glassy systems. Several studies 
have pointed out that the shear stress required to initi- 
ate flow can be larger than the shear stress required to 
prevent slow shear flows from stopping but, does 

this difference in shear stress persist in the infinite system 
size limit? If so, what physical mechanism (for example, 
force chains in granular materials ,16]) is responsible for 
the difference? If not, how significant are the finite-size 
effects? 

We perform molecular dynamics simulations of fric- 
tionless granular materials subjected to boundary-driven 
shear at fixed volume to determine whether or not these 
simple systems can be characterized by a single static 
yield stress in the large system limit. At constant shear- 
ing velocity, we measure the long-time average shear 
stress Tiyv in the limit of zero shearing velocity. We also 
perform simulations at fixed shear force and identify the 
minimum shear stress S^/ required to maintain steady 
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shear flow at long times. We indeed find that E^/ > Ey„ 
at finite system size. However, the difference tends to 
zero in the infinite system size limit. Thus, we argue that 
large frictionless granular systems possess a single static 
yield shear stress. In future studies, we will include static 
friction to determine whether the gap "Eyj — J^yv remains 
finite in large frictional granular systems. 

We also investigate the short-time response of friction- 
less granular systems to applied shear. Previous studies 
of sheared glassy systems ^3 ^ have found that 
the shear stress in response to applied shear strains over- 
shoots the long-time average value at short times. The 
shear stress overshoot is often employed as a dynamic 
measure of the yield shear stress. In addition, these stud- 
ies have found that the size of the overshoot increases 
with increasing shear rate and decreasing temperature. 
Does the shear stress overshoot at short times also occur 
in model granular systems? Is the shear stress overshoot 
related to the difference in the measured values of the 
yield shear stress T,yf — Y,yvl To address these questions, 
we monitor the short-time response of the shear stress to 
applied shear strain over a range of packing fractions and 
shear velocities and in systems where we constrain the ve- 
locity profiles to be linear and in systems without such a 
constraint. We find that the packing fraction and shape 
of the velocity profile strongly influence the short time 
response. In fact, systems near random close packing 
with no constraints on the velocity profile do not possess 
a shear stress overshoot in the range of shear rate consid- 
ered, while systems that are constrained to have linear 
velocity profiles do possess an overshoot. 



II. METHODS 

In this section, we provide important details of the 
simulation methods. We performed molecular dynamics 
simulations of frictionless granular systems in 2D at fixed 
volume in the presence of an applied shear stress. The 
shear stress was applied by moving a top boundary layer 
of particles horizontally as a rigid body at either fixed 
shearing velocity u or fixed lateral force Fq, while the 
bottom boundary remained stationary. We studied sys- 
tems composed of 50-50 mixtures of large and small par- 
ticles with equal mass m and diameter ratio 1.4. These 
bidisperse systems do not crystallize or segregate under 
shear p^.[ig|. 

The position of each particle i in the bulk was ob- 
tained as a function of time t by solving Newton's equa- 
tions of motion 

m-^ = Fi=Y^ [Flj{nj) - b {vi - Vj) ■ hj] fij, (1) 
i 

where the sum over j is a sum over the nearest neighbors 
of particle i. The simple frictionless granular systems 
considered here interact via two pairwise forces that act 
only along the line connecting particle centers fij and are 



nonzero only when particles i and j overlap |20j . The first 
pairwise interaction is the purely repulsive linear spring 
force 

F[^in,)^ — (i-^), (2) 

where e is the characteristic energy scale of the inter- 
action, aij — {ai + crj)/2 is the average diameter of 
particles i and j, and Vij is their separation. The sec- 
ond pairwise interaction is dissipative and proportional 
to velocity differences along f^. We chose the damping 
coefficient b — 0.0375, which corresponds to a restitution 
coefficient e — 0.92 typical for granular systems. 

At constant velocity, the equation of motion for each 
particle in the top boundary is trivial, Cpx/dt^ = 0, sub- 
ject to dx/dt = u. At constant shear force Fq, each par- 
ticle in the top boundary obeys an equation of motion 
similar to that in Eq. ^ 

= ^0 + 5] (F'-(r.j ) ^b{u-v,)- h,) hj ■ i, (3) 

where M is the mass of the top boundary. The second 
term in Eq.|2|is the total horizontal force on particles i in 
the top boundary arising from interactions with particles 
J in the bulk. 

The starting configurations were prepared by choosing 
a packing fraction </> = 0.85 near random close packing 
for this system and random initial particle positions. 
The system was then allowed to relax at fixed volume 
to the nearest local e nerg y minimum using the conju- 
gate gradient method |22|. During the quench, periodic 
boundary conditions were implemented in both the x- 
and 2/-directions. Following the quench, particles with y- 
coordinates y > Ly {y < 0) were chosen to comprise the 
top (bottom) boundary. This preparation algorithm cre- 
ated rough and amorphous top and bottom boundaries, 
which prevents slip between the bulk and boundary par- 
ticles during shear. After the boundaries are constructed, 
the simulation cell was nearly square and contained A^o 
bulk particles and Nt particles in the top and bottom 
boundaries. Periodic boundary conditions in the x direc- 
tion were employed during shear. 

During the simulations, we calculated the shear stress 
on the top and bottom boundaries and in the bulk. Each 
of these measurements gave similar values for the average 
shear stress, however, the shear stress fluctuations were 
much larger on the boundaries as expected. Therefore, 
below we focus on measurements of the bulk shear stress 
calculated using the virial expression p^ : 

^ / JVo ^ No \ 

^ " ~Tj7 { ^ ^'"^■'^'"y' 2 ^ ^v^vv ' (4) 

y \i=l i=^j I 

where bvi — Vi — {vi) is the deviation of the velocity of 
bulk particle i from the average velocity {vi) at height y^. 

We performed several measurements of the shear stress 
near the initiation and cessation of flow. First, at fixed 
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FIG. 1: Deviation in the shear stress E from the shear stress 
Ej,„ in the it — > hmit versus shear rate ujLy for a system 
with A^o ~ 1024 bulk particles sheared at constant u. The 
solid line has slope 0.63. 



shearing velocity, we measured the long-time average 
shear stress Y^y^ in the m ^ limit. Second, at fixed 
lateral force, we measured the minimum shear stress 
S = Fq/Lx required to maintain steady shear flow at 
long times. For all measurements of the shear stress we 
averaged over at least 100 different initial realizations. 
We did not find large differences in the shear stress re- 
sponse among different starting configurations. Also, to 
assess finite size effects, we varied the number of particles 
in the bulk A'o over more than two orders of magnitude 
from Nq — 32 to 4096. In the subsequent discussion of 
results, the small particle diameter a, characteristic en- 
ergy e, and ay/mje were chosen as the units of length, 
energy, and time, and all quantities are normalized by 
these below. 



FIG. 2: Shear strain 7 versus time t for a system with TVo « 
1024 bulk particles sheared at constant force _fo = E / Ly . Two 
shear stresses E = 3.7 x 10"* (solid line) and E = 9.7 x 10"* 
(dashed line) are shown. The smaller value is below and the 
larger is above the minimum yield stress E^/ = 4.4 x 10"* 
required to maintain steady flow at long times. 



form [i[T3 



is the shear stress in the u 



(5) 



where > 0, Y,y^ is the shear stress in the u — > 
limit, and the power-law exponent a k, 0.63 0|. The 
flow curve for a system with TVq « 1024 bulk particles is 
shown in Fig.nand Sj,^ = 2.1 x 10"^ for this system size. 
Systems sheared at constant velocity flow at any nonzero 
u, however, by extrapolating the flow curve to u — > 0, we 
can obtain a measure of the yield shear stress S^^. 



B. Constant Shearing Force 



III. RESULTS 

In this section, we present a number of novel results 
from our simulations of frictionless granular materials 
subjected to boundary-driven planar shear. 



A. Constant Shearing Velocity 

We have measured the average shear stress E as a func- 
tion of shear rate u/ Ly in systems sheared at fixed veloc- 
ity u of the top boundary. At each u, we began with an 
unsheared initial configuration, the system was sheared 
for a strain of at least 10, and then the shear stress was 
averaged over an additional strain of 100. We have shown 
in previous studies that at such large strains these 
systems are spatially uniform and possess linear velocity 
profiles. We find that the flow curve (S versus u/Ly) for 
the system obeys the commonly used phenomenological 



We also studied frictionless granular systems sheared 
at fixed lateral force Fq. In this ensemble, granular sys- 
tems do not flow on long time scales unless the applied 
shear stress E = Fq /L^ exceeds a shear stress threshold. 
In Fig.|21 we show the shear strain 7 — x/ Ly, where 
X is the horizontal displacement of the top boundary, as 
a function of time for applied shear stresses above and 
below Sj,/ in a system with A^o ~ 1024 and averaged 
over 100 initial realizations. When E > S^/, the shear 
strain diverges and the system flows at long times at an 
average shear rate that is consistent with the flow curve 
for the fixed shearing velocity ensemble. The average 
shear rate is given by the slope of strain versus time in 
Fig. 121 When S < Ej,/, the system can flow at short time 
scales. However, the system is able to find a configura- 
tion that can support the applied shear stress, and the 
system stops flowing. Moreover, the flow will not resume 
because dissipation damps the velocity fluctuations. As 
shown in Fig. 13 the maximum shear strain 7s that the 
system attains increases as a power-law with the applied 
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FIG. 3: The difference between the apphed shear stress E and 
the minimum shear stress E.y/ required to maintain shear flow 
at long times versus the maximum shear strain 7^ obtained. 
The sohd line has slope —0.20. 



FIG. 4: The deviation of the yield shear stress Ey from its 
value in the infinite system size limit Eyoo as a function of 
the number of bulk particles A'o in the constant shear velocity 
(circles) and force (squares) ensembles. 



shear stress E and diverges as S ^ Sj,/ from below: 



7. 



,/3' 



(6) 



where Af > 0, Hyf is the minimum shear stress at which 
7s — > 00, and the power-law exponent (3 « 0.20. For 
iVo ~ 1024 bulk particles, S^/ = 4.4 x 10""* > Y.yy. 



Vv 

of S, 



1.1 occurs for iVo > 300. A simple interpretation 



vf 



Tiyy is that 'Syy mcasurcs the average shear 



stress, while Ej,/ is related to the maximum shear stress 
in the u — s- limit. Since the distribution of shear stresses 
becomes a (5-function and the shear fluctuations vanish, 
the difference between T,yy and T,yf vanishes in the large 
system limit for frictionless granular systems. 



C. System-Size Dependence 

In the previous section, we showed that "Eyf > 'Syy for 
systems with A^'o ~ 1024 bulk particles. How does the 
difference in these two measurements of the yield shear 
stress depend on system size? Does the difference tend 
to zero for frictionless granular systems? To answer these 
questions, we performed measurements of Tiyy and Ttyf 
for systems with iVg in the range 32 to 4096. For all 
system sizes, the shear stress obeyed Eq.|Slin the constant 
shear velocity ensemble and Eq. in the constant shear 
force ensemble. We found that T,yf > T,yy for all system 
sizes studied, however, the difference between these two 
measures of the yield shear stress decreased as A'o ^ 00. 
Both measures decreased with increasing system size and 
converged to the same value in the infinite system size 
limit, 'Syoo ~ 1-7 X 10""'. In Fig. ^ we show the system- 
size dependence of "Syf and T,yy. For example, in the 
constant force ensemble, Ej,/ — Ej,oo scales as a power 
law with A'o 

over the entire range of system sizes with Bf > and 
Tjf sa 0.75. Ttyy has a similar power-law dependence for 
small systems, but a more rapid power-law decay with 



D. Discontinuity in the Shear Rate 

In Fig. [SI we show a comparison of the flow curves, 
i.e. shear stress E versus shear rate u/Ly, for systems in 
the constant shearing velocity and force ensembles. For 
large shear stresses E > Ej,/, there is a correspondence 
between shear rate and shear stress in the two ensem- 
bles. However, since E^j is larger than the average shear 
stress T,yy in the m — *■ limit at finite system size, there 
is a jump discontinuity in the shear rate when the ap- 
plied shear stress is increased above E^/. In fact, sev- 
eral recent experimental studies have found that foams, 
emulsions, and granular materials also display a rate of 
strain discontinuity when they begin flowing in response 
to applied shear stress 2&> 22}- However, we find 
that in frictionless granular systems, the jump discon- 
tinuity in shear rate upon initiation of shear flow is a 
flnite size effect — Uc tends to zero in the infinite system 
size limit. The jump discontinuity Uc{No) is obtained by 
solving J^yy{uc, Nq) = T,yf{No). Using the scaling rela- 
tions in Eqs.|31and|7| we find that ^ N^'^''"' - A^g"^'^, 
which is confirmed in the inset to Fig. [S] In light of these 
results, we advocate further experimental studies of the 
rate of strain discontinuities in soft glassy materials, es- 
pecially in planar shear cells, to determine whether there 
are strong finite size effects. 
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FIG. 5: Comparison of the flow curves at constant shear 
velocity (asterisks) and constant shear force (circles) for A^'o ~ 
1024. The solid line is a fit to Eq. |S| Uc/Ly is the jump 
discontinuity in the boundary shear rate that occurs when 
the system begins flowing at constant shear force. The inset 
shows the jump discontinuity Uc/Ly as a function of A'^q. The 
solid line has slope ~ —1.2. 



E. Shear Stress Overshoot 

A frequently used measure of the dynamic yield shear 
stress is the shear stress overshoot above the long-time 
average value in systems sheared at finite shear rate. 
In fact, several recent computational studies of dense 
Lennard- Jones glasses have measured the dependence of 
the shear stress overshoot on the bath temperature and 
imposed shear rate [Tol H IT^ . In this final section, we 
present results from simulations of frictionless granular 
systems undergoing planar shear to determine whether 
a significant shear stress overshoot occurs on short time 
scales in the slowly sheared regime in these systems. 

There are several key differences between our current 
study and previous investigations of the shear stress over- 
shoot in sheared glassy systems: 1) we study systems 
with no constraints on the velocity profile as well as sys- 
tems with constraints that enforce a linear velocity profile 
(vx) = uy/Ly as used in simulations of glasses [23, 2) we 
study a wide range of packing fractions from near random 
close packing at cf) — 0.85 to overcompressed systems at 
(j) = 1.1, and 3) we focus on dissipative, granular systems, 
not conservative, thermal systems. In the results below, 
we show that the packing fraction and shape of the ve- 
locity profile strongly influence the short-time response 
of sheared granular systems. We fixed the dissipation in 
this study, however, the influence of dissipation will be 
investigated in a future study p^ . 

To study the shear stress overshoot, we measured the 
shear stress response of the system at small strains to a 
slow applied shear rate u/Ly = 10~^. The shear stress as 
a function of shear strain averaged over 100 independent 
realizations is shown in Fig. El for a system with iVo w 
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FIG. 6: Shear stress E versus shear strain 7 for systems with 
TVo ^ 4096 at (a) 4> = 0.85 and (b) (j) = 1.10 and u/Ly = 
0.001. The solid and dotted lines show results for systems 
with and without a constraint that enforces linear velocity 
proflles. 



4096. Note that in Fig. O the shear stress is measured 
only over small shear strains up to 7 = 1. In contrast, 
the shear stress was averaged over large shear strains up 
to 7 = 100 in previous figures. 

Two striking results are presented in Fig. First, 
when the velocity profile is unconstrained, the shear 
stress does not overshoot the long-time average shear 
stress at this shear rate. The shear stress increases mono- 
tonically to the long-time average value. The shear strain 
required to reach the long-time average shear stress de- 
creases with increasing 0, but is much less than the shear 
strain required for the system to possess a linear veloc- 
ity profile. See Ref. [13 for an extensive discussion of 
the evolution of the velocity profiles in sheared granular 
systems. In contrast, when we constrain the system to 
possess a linear velocity profile, a large shear stress over- 
shoot develops. In addition, the shear strain required 
to reach the long-time average shear stress is < 0.2 and 
roughly independent of packing fraction. These results 
suggest that the shear stress overshoot at short times in 
glassy and granular systems is an artifact of the fact that 
the velocity field is constrained to be linear. At the very 
least, the constraint significantly amplifies and speeds 
up the response of the system to applied shear. Sec- 
ond, the packing fraction strongly influences the height 
of the shear stress overshoot. The overshoot for (j> = 0.85 
in Fig. (a) (with the constraint) is nearly a factor of 
10 smaller than that for = 1.1 in Fig. (b). More- 
over, if we define A = (E^ — Sp)/Sp as the relative 
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height of the shear stress overshoot with maximum shear 
stress T,jn and shear stress plateau at long times Sp, 
A « 0.09 at = 0.85 compared to A « 0.15 at = 1.10. 
Both of these findings demonstrate that the shear stress 
overshoot is less pronounced in frictionless granular sys- 
tems near random close packing with the local dissipation 
model in Eq. ^ttian in dense Lennard- Jones glasses. 

IV. CONCLUSION 

In this article, we studied model frictionless granular 
systems near the initiation and cessation of shear flow us- 
ing molecular dynamics simulations of boundary-driven 
shear flow at constant volume in 2D. These simulations 
were performed to address several open questions con- 
cerning the jamming (or unjamming) transition in fric- 
tionless granular systems. First, we wanted to determine 
whether these model systems can be characterized by a 
single yield shear stress in the large system limit. We 
compared two measures of yield shear stress: 1) the av- 
erage shear stress Syt, in the limit where the velocity of 
the shearing boundary tends to zero and 2) the mini- 
mum shear stress 'Eyf required for steady shear flow at 
long times when a constant force is applied to the shear- 
ing boundary. As found in previous studies of glassy 
and granular systems, T,yf > Ej,„ in finite-sized sys- 
tems. However, these two measures become identical 
Tiyf = Tiyy 111 thc infinitc system size limit in friction- 
less granular systems. An important direction for future 
research is to determine how the inclusion of static fric- 
tional forces affects these results. Does the difference 
Tjyf — Tiyy TemsLUi finite in the large system limit in fric- 
tional granular systems? Recent studies have argued that 
a large yield stress difference S^/ — E^^ is responsible 
for shear banding — spatially localized velocity profiles — 
in dense Lennard- Jones glasses However, our results 
suggest that another mechanism is responsible for shear 
banding in frictionless granular systems 0, Fur- 
ther studies are required to determine whether a possible 
gap Yiyf — Yiyy > contributes to shear banding in large 



frictional systems. 

Another question addressed in this article is whether 
frictionless granular systems possess a strain rate discon- 
tinuity when they begin flowing at constant force as has 
been found in several recent experimental studies on sim- 
ilar systems [2^ |2^ . We indeed found a discontinu- 
ity in the shear rate upon the initiation of flow, but the 
discontinuity is proportional to T,yf — Y,yy and therefore 
tends to zero in the infinite system size limit. We recom- 
mend further experimental studies in planar shear cells 
to assess the finite size effects on the strain rate discon- 
tinuity. 

We have also investigated the short-time response of 
the shear stress of the system when the shearing bound- 
ary is driven at fixed velocity. It is well-known for glassy 
systems that the shear stress at short times can over- 
shoot the long-time average value when these systems 
are sheared at finite shear rate, and the height of the 
overshoot is often used as a measure of the dynamic yield 
shear stress. We found several novel results for the short- 
time response. First, the shape of the velocity profile 
strongly infiuences the shear stress at short times. When 
our systems did not have a constraint imposed on the 
velocity profile, we did not observe an overshoot in the 
shear stress. In contrast, when a linear velocity profile 
was enforced, a strong shear stress overshoot occurred. 
Second, the height of the overshoot (at least in systems 
with linear velocity profiles) decreases with decreasing 
packing fraction. For example, the height of the peak in 
shear stress is a factor of 10 smaller near random close 
than in an overcompressed system at (f) — 1.1. These re- 
sults point out that there is not a significant shear stress 
overshoot in frictionless granular systems — in contrast to 
dense Lennard-Jones glasses. 
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